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ON THE CONSTRUCTION OF GENERAL CUBATURE 
FORMULA BY FLAT EXTENSIONS 

M. ABRIL BUCERO, C. BAJAJ, AND B. MOURRAIN 


Abstract. We describe a new method to compute general cuba- 
ture formulae. The problem is initially transformed into the com¬ 
putation of truncated Hankel operators with flat extensions. We 
then analyse the algebraic properties associated to flat extensions 
and show how to recover the cubature points and weights from 
the truncated Hankel operator. We next present an algorithm to 
test the flat extension property and to additionally compute the 
decomposition. To generate cubature formulae with a minimal 
number of points, we propose a new relaxation hierarchy of convex 
optimization problems minimizing the nuclear norm of the Han¬ 
kel operators. For a suitably high order of convex relaxation, the 
minimizer of the optimization problem corresponds to a cubature 
formula. Furthermore cubature formulae with a minimal number 
of points are associated to faces of the convex sets. We illustrate 
our method on some examples, and for each we obtain a new min¬ 
imal cubature formula. 


1. Cubature formula 

1.1. Statement of the problem. Consider the integral for a contin- 
nons fnnction /, 

I[f] = [ tn(x)/(x)dx 
Jn 

where C M” and w is a positive fnnction on 

We are looking for a cnbatnre formnla which has the form 

r 

(1) = '^Wj fiCj) 

i=i 

where the points (j G C"^ and the weights Wj G M are independent of 
the fnnction /. They are chosen so that 

W) = i[f],yf ^v, 

where U is a hnite dimensional vector space of fnnctions. Usnally, the 
vector space V is the vector space of polynomials of degree < d, becanse 
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a well-behaved function / can be approximated by a polynomial, so 
that Q[f] approximates the integral /[/]• 

Given a cubature formula Q for I, its algebraic degree is the largest 
degree d for which /[/] = (ct|/) for all / of degree < d. 

1.2. Related works. Prior approaches to the solution of cubature 
problem can be grouped into roughly two classes. One, where the goal 
is to estimate the fewest weighted, aka cubature points possible for sat¬ 
isfying a prescribed cubature rule of hxed degree [911251 [271 EOl [SHIM]. 
The other class, focusses on the determination and construction of 
cubature rules which would yield the fewest cubature points possible 
[3 EH EH SOI SH SH 113 IS]- IHl) example. Radon introduced a 
fundamental technique for constructing minimal cubature rules where 
the cubature points are common zeros of multivariate orthogonal poly¬ 
nomials. This fundamental technique has since been extended by many, 
including for e.g. [Ml EH EH] where notably, the paper [16] use mul¬ 
tivariate ideal theory, while |3^ uses operator dilation theory. In this 
paper, we propose another approach to the second class of cubature 
solutions, namely, constructing a suitable hnite dimensional Hankel 
matrix and extracting the cubature points using sub operators of the 
Hankel matrix [TB]. This approach is related to pil EH Elj and which 
in turn are based on the methods of multivariate truncated moment 
matrices, their positivity and extension properties HHiniiis]. 

Applicatons of such algorithms determining cubature rules and cuba¬ 
ture points over general domains occur in isogeometric modeling and 
hnite element analysis using generalized Barycentric hnite elements 
dH H E3 EZ]. Additional applications abound in numerical integra¬ 
tion for low dimensional (6-100 dimensions) convolution integrals that 
appear naturally in computational molecular biology 13 E], as well in 
truly high dimensional (tens of thousands of dimensions) integrals that 
occur in hnance [33 E]- 

1.3. Reformulation. Let R = ]R[x] be the ring of polynomials in the 
variables x = (xi,..., x„) with coefficients in M. Let Rd be the set of 
polynomials of degree ^ d. The set of linear forms on R, that is, the 
set of linear maps from R to M is denoted by R*. The value of a linear 
form A G R* on a polynomial p G R is denoted by (A|p). The set R* 
can be identihed with the ring of formal power series in new variables 

Z = (Xi, . . .,Zn)-. 

R* M[[z]] 

A ^ A(z) = (A|x“)z“. 
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The coefficients (A|x“) of these series are called the moments of A. The 
evaluation at a point ( G M” is an element of R, denoted by e^, and 
defined by : / e -R ha f{() G M. For any p e R and any A G R*, let 
p-k A q e Re^ A{pq). 

Cubature problem: Let V <Z Rhe the vector space of polynomials 
and consider the linear form J G Id* defined by 

J : 1/ ^ M 

V HA /[v] 

Computing a cubature formula for / on Id then consists in finding a 
linear form 

r r 

a = ■ f 

i=l j=l 

which coincides on Id with I. In other words, given the linear form I 
on Rd, we wish to find a linear form a = Y7i=i which extends I. 

2. Cubature formulae and Hankel operators 

To find such a linear form a & R*, we exploit the properties of its 
associated bilinear form H„ : {p,q) E R x R ^ {o'Ipq), ol equivalently, 
the associated Hankel operator; 

: R -E R* 

p HA px a 

The kernel of is keriLg. = {p G -R | Vg G -R, {cr\pq) = 0}. It is an 
ideal of R. Let Aa = R/kerffa- be the associated quotient ring. 

The matrix of the bilinear form or the Hankel operator associated 
to a in the monomial basis, and its dual are ((A|x"’''^))a^^gNn. If we 
restrict them to a space Id spanned by the monomial basis (x“)Q,g^ for 
some finite set H C N”, we obtain a finite dimensional matrix = 

and which is a Hankel matrix.. More generally, for 
any vector spaces Id, Id' C R, we define the truncated bilinear form 
and Hankel operators: '■ (v,v') G H x Id' ha (cr|pg) G M and 

HJF - Y'*^ jf Y (resp. Id') is spanned by a monomial 

set x"'^ for A C N"'(resp. for B C N"'), the truncated bilinear 
form and truncated Hankel operator are also denoted by The 

associated Hankel matrix in the monomial basis is then [RA^] = 

((A|x“+^))„gA,/3eB- 

The main property that we will use to characterize a cubature for¬ 
mula is the following (see ^31 120] ): 
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Proposition 2.1. A linear form a G R* can be decomposed as a = 

Y7i=i Wi gq with WieC \ {0}, Ci e iff 

• Hfj : p I— )■ p A (T is of rank r, 

• ker is the ideal of polynomials vanishing at the points {Ci, ■ ■ ■ ,Cr 

This shows that in order to find the points Q of a cubature formula, 
it is sufficient to compute the polynomials p & R such that Vg G R, 
{o'\pq) = 0 , and to determine their common zeroes. In section]^ we 
describe a direct way to recover the points Q, and the weights Ui from 
suboperators of H„. 

In the case of cubature formulae with real points and positive weights, 
we already have the following stronger result (see P?1 EU]!: 

Proposition 2.2. Let a E R*. 

r 

a = 

i=l 

with Wi > 0, (i E M” iff rankHf^ = r and 0. 

A linear form a = Wi > 0, G M"’ is called a r- 

atomic measure since it coincides with the weighted sum of the r Dirac 
measures at the points Q. 

Therefore, the problem of constructing a cubature formula a for I 
exact on V <Z R can be reformulated as follows: Construct a linear 
form a E R* such that 

• ranki^CT = r < oo and ^ 0. 

• \/v eV, I[v] = (cT|n). 

The rank r of is given by the number of points of the cubature 
formula, which is expected to be small or even minimal. 

The following result states that a cubature formula with dim(I/) 
points, always exist. 

Theorem 2.3. [12111] If a sequence (cTa)Q,gr!}n is the truncated mo¬ 
ment sequence of a measure mu (i.e. o'a = f x°‘dmu for |a:| ^ t) 
then it can also be represented by an r-atomic measure: for |q;| ^ t, 

= Y7i=iWi^f where r ^ St, Wi > 0, Q E supp{mu). 

This result can be generalized to any set of linearly independent 
polynomials vi,...,Vr G R (see the proof in |1] or Theorem 5.9 in 
1231). We deduce that the cubature problem always has a solution with 
dim(I/) or less points. 

Definition 2.4. Let rc{I) be the maximum rank of the bilinear form 
hY'^' : {w,w') eW xW ^ I[ww'] where W, W' C P are such that 
\fw E W,\/w' E W, ww' eV . It is called the Catalecticant rank of I. 
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Proposition 2.5. Any cuhature formula for I exact on V involves at 
least rc{I) points. 


Proof. Suppose that cr is a cubature formula for I exact on V with r 
points. Let W, W C P be vector spaces are such that Vw G hh, Vtc' G 
hh', ww' . Since hY’^ coincides with HY’^' , which is the restric¬ 
tion of the bilinear form Her to W x W', we deduce that r = rankifer ^ 
rankiff'^' = rank(i7f’'^'). Thus r ^ re(I). □ 


Corollary 2.6. Let IV C V such that \/w,w' G W,ww' G V. Then 
any cubature formula of I exact on V involves at least dim(iy) points. 


Proof. As we have Vp G W, G C so that /(p^) = 0 implies p = 0. 


Therefore the quadratic form jjY’^ : {p, q) E W xW ^ ^\pq] is positive 
definite of rank dim(lT). By Proposition 2.5| a cubature formula of I 
exact on V involves at least rc{I) ^ dim(lT) points. □ 


In particular, if 1/ = Rd any cubature formula of / exact on V involve 
at least dimi?|^dj = points. 

In PS], this lower bound is improved for cubature problems in two 
variables. 


3. Flat extensions 

In order to reduce the extension problem to a finite-dimensional prob¬ 
lem, we consider hereafter only truncated Hankel operators. Given 
two subspaces W, W of R and a linear form a defined on W ■ W (i.e. 
a G {W ■ IT')*), we define 

HY’^' : it X it' ^ M 

{w,w') H-)■ {a\ww'). 

If w (resp. w') is a basis of IT (resp. IT'), then we will also denote 
^w,w' ._ ffW,w\ matrix of in the basis w = {wi, ..., Wg}, w' = 

{w'l,..., is 

Definition 3.1. Let IT C T, IT' C V' be subvector spaces of R and a G 
{V-V'Y. ITe say that HY^' is a flat extension of HY’^' if rankHY’^' = 
ranki/f’^'. 

A set B of monomials of R is connected to 1 if it contains 1 and if 
for any m Y ^ ^ Bi there exists 1 < i < n and m' E B such that 
m = Xim'. 

As a quotient i?/ker Q has a monomial basis connected to 1, so in 
the first step we take for w, w', monomial sets that are connected to 1 
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For a set B of monomials in i?, let us define = BVJxiBU- ■ -UXn, B 
and dB = B+\ B. 

The next theorem gives a characterization of flat extensions for Han- 
kel operators dehned on monomial sets connected to 1. It is a general¬ 
ized form of the Curto-Fialkow theorem |13) . 

Theorem 3.2. Let B a C, B' G C be sets of monomials 

connected to 1 such that \B\ = \B'\ = r and C ■ C contains B^ ■ B'~^. 
If cr G {C ■ C'Y is such that rankif^’^' = rankif^’*^' = r, then H^’^' 
has a unique flat extension for some a E R*. Moreover, we have 
kerif^ = and R = {B) ©kerif^ = {B') ©kerifg.. In the 

case where B' = B, if 0, then 0. 

Based on this theorem, in order to hnd a flat extension of , it 

suffices to construct an extension of the same rank r. 

Corollary 3.3. Let V G R be a finite dimensional vector space. If 
there exist a set B of monomials connected to 1 such that V G {B'^-B^) 
and a G {B~^ ■ B^)* such that Vu G V, (ulu) = I[v] and rankhf^’^ = 
lankllY^'^^ = \B\ = r, then there exists Wi > 0, Ci ^ such that 
Wv G V, 

r 

i=l 

This characterization leads to equations which are at most of de¬ 
gree 2 in a set of variables related to unknown moments and relation 
coefficients as described by the following proposition; 


Proposition 3.4. Let B and B' be two sets of monomials of Rof sizer, 
connected to 1 and a be a linear form on {B'^ ■ B'^). Then, 
admits a flat extension such that is of rank r and B (resp. B’) 
a basis of R/ kei iff 


( 2 ) 


[H^ 


B+,B'+-i _ 



M' 

N 




with Q = M' = M* = N = is such 

that Q is invertible and 

(3) M = Q*P,M'= QP',N = P*QP', 

for some matrices P G , P' G 


Proof. If we have M = Q* P, M' = Q P', N = P* Q P', then 
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3.2 


has clearly the same rank as Q=[H^’^']. According to Theorem 
^B+,s'+ flgi; cxtenslon Ha with a E R* snch that B and B' 


are bases of = R/ ker ifg.. 

Conversely, if is a flat extension of with B and B' bases 

of Aa = R/ kei Ha, then [H?’^ ] = [H^’^'] = Q is invertible and of size 
r = \B\ = \B'\. As Ha is of rank r, is also of rank r. Thns, 

there exists P' G (P' = Q“^M') snch that M' = QP'. Similarly, 

there exists P G snch that M = Q*P. Thns, the kernel of 

[^iJ+,B'+] (resp, = is the image of ^ ^ (resp. 

^)- We dednce that N = M^P' = P^QP'. □ 


Remark 3.5. A basis of the kernel of is given by the columns 

f -P' \ 

of { j 1, which represent polynomials of the form 


Pa = X“ - ^ Pa,l3 

y&B 

for a G dB. These polynomials are border relations which project the 
monomials x" of dB on the vector space spanned by the monomials 
B, modulo ker It is proved in [B] that they form a border 

basis of the ideal ker Ha when is a flat extension and Hf’^' is 

invertible. 


Remark 3.6. Let A C N"" be a set of monomials such that {cr|x") = 
/[x"]. Considering the entries o/P, P' and the entries a a of Q with 
a ^ A as variables, the constraints are multilinear eguations in 
these variables of total degree at most 3 if Q contains unknown entries 
and 2 otherwise. 


Example 3.7. fTe consider here V = i? 2 fc for k > 0. By Proposition 
3 . 4 , any cubature formula for I exact on V has at least r^ := dimR^. 


points. Let us take B to be all the monomials of degree < k so that B^ 
is the set of monomials of degree < k + 1. If a cubature formula for I 
is exact on i? 2 fc and has r^ points, then is a flat extension of 

Hf'^ of rank r^. Consider a decomposition of as in §). By 

we have the relations 


Proposition jS.f 

(4) 


= Q P, N = P* 


where 
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• M = ((cr|x^+^'))/ 3 gs_/ 3 /gaB with (ct|x^+^') = /[x^+^'] when \(5 + 

13'\ < 2 k, 

• N = ((a|x^+^'))^,/3/6as 

• P = {pi3,a)l3eB,aedB- 


The equations 0 are quadratic in the variables P and linear in the 
variables in M. Solving these equations yields a flat extension 
of . As 0, any real solution of this system of equations 

corresponds to a cubature for I on exact i? 2 fc of the form a = ^0 

with Wi > 0, Ci ^ 

We illustrate the approach with R = M[a:i, X 2 \,V = R^, B = {1, xi, X 2 ,x\, XiX 2 ,xf}, 
= {l,Xi,X 2 ,x\,XiX 2 ,x\,x\,x\x 2 ,Xix\,x\}. Let 


a = 8 + 17 22 — 4 zi + 15 ^2 + 14 ziZ2 — W z\ + 431 z\ — Q + 34 zlz2 — 52 zf 
+51 4 + 38 zizi - 18 zfzi + 86 zfz2 - 160 zf 


be the series truncated in degree 4, corresponding to the first moments 
(not necessarily given by an integral). 


H. 


B+,B+ _ 


8 

-4 

17 

-16 

14 

15 

-52 

34 

-6 

47 

-4 - 

-16 

14 

-52 

34 

-6 

-160 

86 

-18 

38 

17 

14 

15 

34 

-6 

47 

86 

-18 

38 

51 

-16 - 

-52 

34 

-160 

86 

-18 

(7l 

<72 

<73 

(74 

14 

34 

-6 

86 

-18 

38 

(72 

<73 

(74 

<75 

15 

-6 

47 

-18 

38 

51 

(73 

(74 

(75 

<76 

-52 - 

-160 

86 

C7l 

(72 

£73 

<77 

<78 

(79 

(7l0 

34 

86 

-18 

£72 

(73 

(74 

(78 

<79 

<710 

(711 

-6 - 

-18 

38 

(73 

(74 

£75 

(79 

(7l0 

(711 

(712 

47 

38 

51 

(74 

(75 

£76 

<710 

(711 

<712 

(713 

<75,05 <72 

= 

<74,1, (72 

= <73,2 

5 <74 

= <72,3 

5(75 = 

<7i,45 <76 


<7o,5? 

= <75,1 

5 <79 

= <74,25 

<7l0 = 

CO 

CO 

<7ll = 

<72,45 <712 — <7l 

,55 (7i3 — 


< 70,6 ■ 

The first 6x6 diagonal block is invertible. To have a flat 

extension , we impose the condition that the sixteen 7x7 minors 

of , which contains the first 6 rows and columns, must vanish. 

This yields the following system of quadratic equations: 
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- -814592 a1 - 1351680 aicra - 476864 aias - 599040 ai - 301440 aacrs - 35072 cr| 

-520892032(71 - 396821760(72 - 164529152 as + 1693440(77 - 86394672128 = 0 
-814592 ai - 1351680 a2(73 - 476864(72(74 - 599040(71 - 301440 a3(74 - 35 0 72 (7| 

+335275392 (72 + 257276160 (73 + 96277632 (74 + 1693440 (79 - 34904464128 = 0 
-814592(71 - 1351680 a3(74 - 47 6 864(73(75 - 59 9 0 40 (7| - 301440 a4(75 - 35 0 72 erf 
+13226880(73 + 13282560(74 - 8227200(75 + 1693440 an + 31714560 = 0 
-814592 a| - 1351680 a 4 a 5 - 476864 aiag - 599040 af - 301440 asag - 35072 a| 

+212860736 ai + 162698880 a5 + 51427456 ag + 1693440 ais - 13356881792 = 0 
-814592 aia2 - 675840 aia3 - 238432 aia4 - 675840 a^ 

-837472 a 2 a 3 - 150720 a 2 a 4 - 150720 a| - 35072 a 3 a 4 + 167637696 ai - 131807936 aa 
-150272064 a3 - 82264576 a4 + 1693440 ag + 54990043008 = 0 
-814592 aaag - 675840 a 2 a 4 - 238432 a 2 a 5 - 675840 a| - 837472 a 3 a 4 - 150720 a 3 a 5 - 150720 a| 

-35072 a4a5 + 6613440 aa + 174278976 aa + 124524480 a4 + 48138816 aa + 1693440 aio - 746438400 = 0 
-814592 a 3 a 4 - 675840 agas - 238432 aaag - 675840 a| - 837472 a 4 a 5 - 150720 a 4 a 6 - 150720 a| 

-35072 aaae + 106430368 aa + 87962880 a4 + 32355008 as - 4113600 ag + 1693440 aia - 183201600 = 0 
-814592 a 2 a 4 - 675840 aaas - 238432 aaag - 675840 asaa - 599040 agas - 150720 agag - 238432 a| 

-150720 a 4 a 5 - 35072 a 4 a 6 + 106430368 ag + 81349440 ag + 193351424 ag + 128638080 ag + 48138816 ag 
+ 1693440 an - 21721986624 = 0 

-814592 aiag - 675840 ajag - 238432 ajag - 675840 agag - 599040 agag - 150720 agag - 238432 a| 
-150720 agag - 35072 agag + 6613440 ag + 6641280 aj - 264559616 ag - 198410880 ag - 82264576 ag 
+ 1693440 ag + 1312368000 = 0 

—814592 agag — 675840 agag — 238432 agag — 675840 agag — 599040 agag — 150720 agag — 238432 agag 
-150720 agag - 35072 agag + 106430368 ag + 81349440 aa + 25713728 ag - 260446016 ag 
-198410880 ag - 82264576 ag + 1693440 ago + 34550702464 = 0 


The set of solutions of this system is an algebraic variety of dimen¬ 
sion 3 and degree 52. A solution is 

ag = —484, aa = 226, ag = —54, ag = 82, ag = —6, ag = 167, ar = —1456, ag = 614, 
ag = —162, ago = 182, agg = —18, aga = 134, agg = 195. 

3.1. Computing an orthogonal basis of Aa-- In this section, we 
describe a new method to construct a basis B of and to detect flat 
extensions, from the knowledge of the moments (Tq of ct(z). We are 
going to inductively construct a family P of polynomials, orthogonal 
for the inner product 


(p,g) ^ {p,q)a := {a \ pq), 

and a monomial set B connected to 1 such that {B) = (P). 

We start with B = {1}, P = {1} C R. As (1, l)^- = ((+ | 1) 7 ^ 0, the 
family P is orthogonal for a and (B) = (P). 

We now describe the induction step. Assume that we have a set 
B = {mi,... ,ms} and P = [pi,... ,ps} such that 

• (B) = {P) 

• 7 ^ 0 if i = j and 0 otherwise. 
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To construct the next orthogonal polynomials, we consider the mono¬ 
mials in dB = {m'l ,..., and project them on (P); 


P'i = K 


S 

j 


,rn. 


i,Pj)c 


{Pj:Pj)c 


-Pj,i = 


By construction, {p'i,Pj)a = 0 and {pi, ■ ■ ■ ,Ps,p'i) = (mi, ..., m'). 
We extend B by choosing a subset of monomials B' = {m'^, • • •, m'^} 
such that the matrix 

[{Pij ’ Pij/ )O'] 

is invertible. The family P is then extended by adding an orthogonal 
family of polynomials {p^+i,... ,Ps+k} constructed from {p'^,... 

If all the polynomials p' are such that (p(,p()o- = 0, the process stops. 
This leads to the following algorithm; 


Algorithm 1. 

Input: the coefficients a a of a series a G C[[z]] for a G A C N” 
connected to 1 with ao ^ 0. 

• Let B := {!}; P = {!}; r .•= 1; E = (z")aeA; 

• While s > 0 and B^ ■ P+ C E do 

- Compute <9P = {m'l,..., m(} andp' = 

— Compute a (maximal) subset B' = {m'^,..., m'^} of dB 
such that [{Pi.,p'i.,)a-]i!^j,j%k is invertible. 

— Compute an orthogonal family of polynomials {p^+i,... ,Ps+k 
from {p'^,...,p'J. 

- P ;= P U P', P := P U {p^+i,... ,Ps+k}; r+= k; 

• If B'^ ■ B^ (fi E then return failed. 

Output: failed or success with 

• a set of monomials B = {mi,..., m,.} connected to 1, and non¬ 
degenerate for (•,-)cr. 

• a set of polynomials P = (pi ,... ,pr} orthogonal for a and such 
that (P) = (P). 

• the relations pi := m' — monomials m' 

in (9P = (m'l,..., m(}. 

The above algorithm is a Gramm-Schmidt-type orthogonalization 
method, where, at each step, new monomials are taken in dB and 
projected onto the space spanned by the previous monomial set P. 
Notice that if the polynomials pi are of degree at most d' < d, then only 
the moments of a of degree ^ 2 d' -|- 1 are involved in this computation. 
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Proposition 3.8. If Algorithm outputs with success a set B = 
{mi,... ,mr} and the relations pi := m' — Pj> 

dB = (m'l,... ,rn[}, then a coincides on {B^ ■ B^) with the series a 
such that 


• ranki^CT = t ; 

• B and P are basis of Aa for the inner product (•, ■)„; 

• The ideal la = kerif^ is generated by 

• The matrix of multiplication by Xk in the basis P of Aa is 


Mu 


{a I XkPiPj) \ 


Proof. By construction, B is connected to 1. A basis B' of {B~^) is 
formed by the elements of B and the polynomials pi,i = 1,... ,1. Since 
Algorithmic stops with success, we have Vz,j G [1,^], V6 G {B), 
{phb)a = {pi:pj)a = 0 aud pi,...,p; G As (5+) = 

(B) © (pi,... ,p;), = rankif®’^ and is a flat ex¬ 

tension of . By construction, P is an orthogonal basis of {B) and 
the matrix of in this basis is diagonal with non-zero entries on 

the diagonal. Thus is of rank r. 

By Theorem 3.2, a coincides on {B^ ■ B^) with a series d G i?* such 
that i? is a basis of Aa = R/Ia and la = (kerhf^ ’ ) = (pi ,... ,pi). 

As (5+) = (B) © (pi ,.. .,pi) = (P) © (pi ,... ,pi) and P is an or¬ 
thogonal basis of Aa, which is orthogonal to (pi,..., p;), we have 


^kPi = 

t=i 


(a I XkPiPj) 

I P]) 


Pj + P 


with p G (pi,... ,pi). This shows that the matrix of the multiplication 
by Xu modulo la = (pi,..., pi), in the basis P = (pi,... ,Pr} is Mu = 


( A\^kPiPj) \ 

\~Rpfr 


□ 


Remark 3.9. It can be shown that the polynomials (pj)j=i,...,i are a 
border basis of la for the basis B jMl 111 1211 [29] . 


Remark 3.10. If 0, then by Proposition 2.2, the common 

roots Qi,... Ar of the polynomials pi,..., p/ are simple and real G MP. 
They are the cubature points: 


i=l 


with Wj > 0. 
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4. The cubature formula from the moment matrix 


We now describe how to recover the cubature formula, from the 
moment matrix . We assume that the flat extension condition 

is satisfied: 


(5) 


rank H, 




= = \B\ = \B' 


Theorem 4.1. Let B and B' he monomial subsets of R of size r 
connected to 1 and a G {B~^ ■ B'~^)*. Suppose that = 

= \B\ = |5'|. Let Mi = and (M/)* = 

Then, 

(1) B and B' are bases of Aa = R/kerH^, 

(2) Mi (resp. M-) is the matrix of multiplication by x, in the basis 
B (resp. B') of A^, 


Proof. By the flat extension theorem 3^, there exists a G R* such that 
is a flat extension of of rank r = \B\ = \B'\ and kerif^. = 

(ker As R = (B) ©kerif^- and rankhfs- = r, A^ = R/kerH^ 

is of dimension r and generated by B. Thus B is a basis of Aa. A 
similar argument shows that B' is also a basis of Aa. We denote by 
TT : Aa —t (B) and vr' : Aa — t (B') the isomorphism associated to these 
basis representations. 

The matrix ] is the matrix of the Hankel operator 

Ha : Aa —t A^a 
a I—)■ a A ct 


in the basis B and the dual basis of B'. Similarly, is the 

matrix of 

H -■ A- — A* 

a HA a-k Xi'kd 

in the same bases. As Xj a ct = do Mi where Mi : Aa —t Aa is 
the multiplication by Xj in Aa, we deduce that Hx^^a = Ha o Mi and 
] jg matrix of multiplication by Xj in the basis 
B of Aa- By exchanging the role of B and B' and by transposition 
{[H^’^'Y = [H^''^]), we obtain that ^^e transpose 

of the matrix of multiplication by Xj in the basis B' of Aa- D 

Theorem 4.2. Let B be a monomial subset of R of size r connected to 
1 and a G {B^-B^)*. Suppose that lankH^^’^^ = rankif^’-® = |i?|=r 
and that H^’^ ^ 0. Let Mi = Then a can he 
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decomposed as 

r 

i=i 

with Wj > 0 and Q G M” such that Mi have r common linearly inde¬ 
pendent eigenvectors Uj, j = 1,..., r and 

• = 7 ^^ forl<i<n,l<j <r. 


w = __. 


Proof. By theorem 4^, the matrix Mi is the matrix of multiplication 
by Xi in the basis B of A^. As 0, the flat extension theorem 


3.2 implies that 0 and that 


i=i 


where Wj > 0 and Q G M"" are the simple roots of the ideal kerif^- 
Thus the commuting operators Mi are diagonalizable in a common 
basis of eigenvectors Uj, i = 1,..., r, which are scalar multiples of the 
interpolation polynomials at the roots Ci) • • • Xr'- Ui(Ci) = 7 ^ 0 and 

Uj(Cj) = 0 if j 7 ^ i (see [T5j[Chap. 4] or |ini)- We deduce that 

r 

{a I Uj) = '^WkUj{Ck) = WjXj and (d | XiUj) = Q^iWjXj, 

k=l 


so that (j^i 


(g-diUj) 

(o-|u/) 


As Ui(0) 


Aj, we have Wj 


(g'luj) 


□ 


Algorithm 2. 

Input: B is a set of monomials connected to 1, cr G {B~^ ■ B^)* such 
that is a flat extension of of rank \B\. 

• Compute an orthogonal basis {pi,... ,pr} of B for a; 


Compute the matrices M^ = ^ j 




• Compute their common eigenvectors u^,..., u^; 
Output: For j = 1,... ,r, 


Q = 


Wi = 


(ctIxiUj) {a\xinVLj) 

, . . . , 


(o-|Uj) 

Huj) 


(fT|Uj> 




Remark 4.3. Since the matrices commute and are diagonalizable 
with the same basis, their common eigenvectors can be obtained by com¬ 
puting the eigenvectors of a generic linear combination liMi + ■ ■ ■ + 
InMn, k G M. 
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5. CUBATURE FORMULA BY CONVEX OPTIMIZATION 

As described in the previous section, the computation of cubature 
formulae reduces to a low rank Hankel matrix completion problem, 
using the flat extension property. In this section, we describe a new 
approach which relaxes this problem into a convex optimization prob¬ 
lem. 

Let V C R he a vector space spanned by monomials x" for a G A C 
Our aim is to construct a cubature formula for an integral function 
I exact on V. Let i = (/[x"])Qg^ be the sequence of moments given by 
the integral I. We also denote i G Id* the associated linear form such 
that Vn e Id{i I n) = I[v]. 

For A; G N, we denote by 

■^*^(1) = {FTcr I O' G Rif,, (Ja = Iq for a ^ A, 0}, 

the set of semi-dehnite Hankel operators on Rt is associated to moment 
sequences which extend i. We can easily check that 'H^(i) is a convex 
set. We denote by Rri^) the set of elements of of rank ^ r. 

A subset of 'H^(i) is the set of Hankel operators associated to cuba¬ 
ture formulae of r points: 

= I iL. G n^i) I a = > 0,0 G M" 

[ i=i 

We can check that (i) is also a convex set. 

To impose the cubature points to be in a semialgebraic set S defined 
by equality and inequalities 5 = {x G M"' | 5 '?(x) = 0,... = 

0 ,( 7 j*“(x) ^ 0,..., ( 702 (^) ^ 0}; can rehne the space of by 

imposing that a is positive on the quadratic module (resp. preordering) 
associated to the constraints g = {gi,..., g^', gi ,..., jl9]. For the 

sake of simplicity, we don’t analyze this case here, which can be done 
in a similar way. 

The Hankel operator H„ G T^(i) associated to a cubature formula 
of r points is an element of 'H^(i). In order to hnd a cubature formula 
of minimal rank, we would like to compute a minimizer solution of the 
following optimization problem; 

min rank(iL) 

However this problem is NP-hard |16| . We therefore relax it into the 
minimization of the nuclear norm of the Hankel operators, i.e. the 
minimization of the sum of the singular values of the Hankel matrix 
|38| . More precisely, for a generic matrix P G we consider the 
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following minimization problem; 

(6) min trace(P*iJP) 

Let {A, B) E X —)■ {A, B) = tTace{AB) denote the inner 

product induced by the trace on the space of Sk x Sk matrices. The 
optimization problem (1^ requires minimizing the linear form H —)■ 
trace(PPP^) = {H,PP^) on the convex set 'P^(i). As the trace of 
P^HP is bounded by below by 0 when if 0 , our optimization 
problem ([^ has a non-negative minimum ^ 0. 

Problem (|^ is a Semi-Dehnite Program (SDP), which can be solved 
efficiently by interior point methods. See |S|. SDP is an important 
ingredient of relaxation techniques in polynomial optimization. See 
[191123]. 

Let = |p = I Pi € Pfcj be the set of polynomials of de¬ 

gree ^ 2 k which are sums of squares, let be the vector of all 
monomials in x of degree ^ k and let g(x) = (x*^^^)*PP* x*^^) E 

Let Pi(x) (1 ^ z ^ Sk) denote the polynomial (Pj,x*^^)) associated to 
the column Pj of P. We have g(x) = ^ ^ 

trace(P*P,P) = {H, \ PP*) = {a \ g(x)). 

For any / G M, we denote by tt; : P; —)■ Ri the linear map which 
associates to a polynomial p E Ri its homogeneous component of degree 
1. We say that P is a proper matrix if 7r2fc(q'(x)) ^ 0 for all x G 

We are thus looking for cubature formula with a small number of 
points, which correspond to Hankel operators with small rank. The 
next results describes the structure of truncated Hankel operators, 
when the degree of truncation is high enough, compared to the rank. 

Theorem 5.1. Let a E he its truncated Hankel operator 

on Rk- If H(j is of rank r ^ k, then 

r' T 

cr = ^Wje^, ^ Wie^, o 712^ 

i=l i=r'-|-l 

with Wi G C \ {0}, Ci G C” distinct. If moreover 0, then oji > 0 

and (i E MA for i = 1,... ,r. 

Proof. The substitution tq : S'pfc] —)• R 2 k which replaces xq by 1 is an 
isomorphism of K-vector spaces. Let Tq : P^^ —)■ be the pull-back 
map on the dual (TQ(cr) = a o tq). Let a = a o tq = Tq (a) G S'^^k] 
the linear form induced by a on S'[ 2 fc] and let : S[k] —t be the 

corresponding truncated operator on S^k]- The kernel K of is the 
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vector space spanned by the homogenization in xq of the elements of 
the kernel K of Hu. 

Let be the lexicographic ordering snch that xq >“ • ■ ■ >“ Xn- By 
114| [Theorem 15.20, p. 351], after a generic change of coordinates, the 
initial J of the homogeneons ideal {K) C S' is Borel fixed. That is, if 
XjX" G J, then xjx" G J for j > i. Let B be the set of monomials 
of degree k, which are not in J. Note J is Borel fixed and different 
from S'[ 2 fc], Xq G B. Similarly we check that if Xq° 


ai-\ 


X 




G J. 


x;;" G J with 
This shows that 


ai = ■ ■ ■ = ai-i = 0, then Xg 
B = tq{B) is connected to 1. 

As {B) © (J) = {B) ® K = S'[fc] where K = ker we have \B\ = r. 
As B is connected to 1, deg(i?) < r ^ k and B~^ C Rk- 

By the snbstitntion xq = 1, we have Rk = (B) (B K with K = 
ker Hu- Ther efore, is a flat extension of H^’^. Via the flat extension 
theorem 3.2, there exists A* G IK \ {0}, Q = ..., G IK"'”''^ 

distinct for z = 1,..., r snch that 


(7) 


a = J2 on Si 


[2k]- 


2 = 1 


Notice that for any A 7 ^ 0, on S'[ 2 fc]. 

By an inverse change of coordinates, the points Q of Q are trans¬ 
formed into some points 0 = ■ ■ ■ Xi,n) ^ IK"'” such that 

0,0 7 ^ 0 (say for i = 1 ,... ,r') and the remaining r — r' points with 
Cj,o = 0. The image by Tq of e^. G S'[ 2 fc] with 7 ^ 0 is 

^■0(00) = 

where 0 = ^(0,i, • • •, 0,n)- The image by of G S'i^ 2 fc] ^i^h 0,o = 0 
is a linear form G i? 2 fc 5 which vanishes on all the monomials x“ with 
|a| < 2 k, since their homogenization in degree 2 k is Xq*^~'“'x“ and their 
evaluation at Q = (0, Ci,i, • • •, 0,n) gives 0. The value of TQ(e^J at 
x“ with |a| = 2k Q ■ ■ ■ where 0 = (0,i, • ■ ■, Ci,nj- We 

deduce that 

T'o(ed) = 

The flat extension theorem implies that if H^ ^ 0 then Aj > 0 and 
(i G in the decomposition Q. By dehomogenization, we have 

> 0, 0 = 7^(0,1, • • • ,0,n) e M" for z = l,...,r' and 

□ 


0 = (0,1, • • •, 0,n) G M" for z = r' + 1,..., n. 


We exploit this structure theorem to show that if the truncation 
order is sufficiently high, a minimizer of (|^ corresponds to a cubature 
formula. 
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Theorem 5.2. Let P be a proper operator and k ^ Assume 

that there exists a* G be such that H„* is a minimizer of (0 of 

rank r with r ^ k. Then H^* G Sf{i) i.e. there exists cuj > 0 and 
(fi G M" such that 

r 

a* = 

i=l 


Proof. By Theorem 5.1 


r' r 

O'* = ^ O TTsfc, 

i=l i=r'+l 


with (Ui > 0 and G M"" for i = 1,..., r. 

Let us suppose that r 7 ^ r'. As fc ^ deg(u)+i ^ elements of V are 

of degree < 2k, therefore a* and a' = coincide on V and 

H^i G We have the decomposition 


trace(Pi/<^*P) = {a* \ q) = {a \ q)+'^ 0 JiT^ 2 k{q){Ci)- 

i=l 

The homogeneous component of highest degree 7 i2k{q) of g(x) = Piii^Y 
is the sum of the squares of the degree-fc components of the pp. 

Sk 

'^2k{q) = 5^(7rfc(Pi))^ 
i=l 

so that 'Yl\=i 00 iT^ 2 k{q){Ci) ^ 0- As trace(PPo-*P) is minimal, we must 
have YJi=r'+i 0 ^i'^ 2 k{,q){,C,i) = 0, which implies that 7 r 2 A:(g)( 0 ) for i = 
r' + 1,... ,r. However, this is impossible, since P is proper. We thus 
deduce that r' = r, which concludes the proof of the theorem. □ 


This theorem shows that an optimal solution of the minimization 
problem ([^ of small rank (r ^ k) yields a cubature formula, which is 
exact on V. Among such minimizers, we have those of minimal rank 
as shown in the next proposition. 

Proposition 5.3. Let k ^ and H be an element of'H'^(i) with 

minimal rank r. If k ^ r, then H G Srif) ond it is either an extremal 
point o/P^(i) or on a face of'H^{i), which is included in T^(i). 


Proof. Let Po- € £. 


By Theorem (5.1 
and (i G ^ 


be of minimal rank r. 


) O' = Y.l=i + El=r '+1 ° with Ui>0 


for i = 1,..., r. The elements of V are of degree < 2k, 
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therefore a and a' = Y7i=i coincide on V. We deduce that 
Ha' e •H'^(i). 

As ranki/o-/ = r' ^ r and H„ G 'H^(i) is of minimal rank r, r = r' 
and G 

Let us assume that is not an extremal point of 'H^(i). Then it is 
in the relative interior of a face F of 7/^(1). For any in a sufficiently 

small ball of F around H^, there exist t g] 0, 1[ and G F such that 

Ha = tHai + (1 — t)Ha2- 

The kernel of H^ is the set of polynomials p E Rk such that 

0 = Ha{p,p) = tHa,{p,p) + {1 -t)Ha2{p,p). 

As iLo-i ^ 0, we have Ha^{p,p) = 0 for i = 1,2. This implies that 
keriLo- C keriLo-i, for * = 1;2. From the inclusion keriLo-^ PikeriLo -2 C 
ker Ha, we deduce that 

ker Ha = ker Ha^ P ker Ha^ ■ 

As Ha is of minimal rank r, we have dimkeriLg. ^ dimkeriLo-r This 
implies that kei Ha = kei Ha^ = keiHa 2 - 

As r ^ k, Aa has a monomial basis B (connected to 1) in degree 
< k and Rk = {B) © kerFfo-- Consequently, Ha (resp. Ha^) is a flat 
extension of H^’^ (resp. H^.’^) and we have the decomposition 

r r 

i=l i=l 

with cjjj > 0, i = 1,..., r, j = 1, 2. We deduce that Ha^ £ 

all the elements of the line {Ha, Ha^) which are in F are also in £r{^)- 

Since F is convex, we deduce that F <Z £^{\). □ 

Remark 5.4. A cuhature formula is interpolatory when the weights are 
uniquely determined from the points. From the the previous Theorem 
and Proposition, we see that if a cubature formula is of minimal rank 
and interpolatory, then it is an extremal point of'H^{i). 

According to the previous proposition, by minimizing the nuclear 
norm of a random matrix, we expect to hnd an element of minimal 
rank in one of the faces of HA{\), provided k is big enough. This yields 
the following simple algorithm, which solves the SDP problem, and 
checks the flat extension property using Algorithm Furthermore it 
computes the decomposition using Algorithm or increases the degree 
if there is no flat extension: 


Algorithm 3. 
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k : = 


deg(U) 

2 


; notflat := true; P:= random Sk x Sk matrix; 

While (notflat) do 

— Let a be a solution of the SDP problem: trace(P*iJP); 

— If is not a flat extension, then k := k + 1; else not- 
flat:=false; 

Compute the decomposition of a = >0,0^ 


6. Examples 

We now illustrate our cubature method on a few explicit examples. 


Example 6.1 (Cubature on a square). Our first application is a well 
known case, namely, the square domain 12 = [—1,1] x [—1,1], We solve 
the SDP problem with a random matrix P and with no constraint 
on the support of the points. In following table, we give the degree 
of the cubature formal (i.e. the degree of the polynomials for which 
the cubature formula is exact), the number N of cubature points, the 
coordinates of the cubature points and the associated weights. 


Degree 

N 

Points 

Weights 

3 

4 

±(0.46503, 0.464462) 
±(0.855875, -0.855943) 

1.545 

0.454996 

5 

7 

±(0.673625, 0.692362) 
±( 0 . 40546 , -0.878538) 
±(-0.901706, O. 34 O 6 I 8 ) 
(0, 0) 

0.595115 

0.43343 

0.3993 

1.14305 

7 

12 

±(0.757951, 0.778815) 
±(0.902107, 0.0795967) 

±(0.04182, 0 . 9432 ) 

±(0.36885, 0 . 19394 ) 
±(0.875533, -0.873448) 
±(0.589325, - 0 . 54688 ) 

0.304141 

0.203806 

0.194607 

0.756312 

0.0363 

0.50478 


The cubature points are symmetric with respect to the origin (0, 0). 
The computed cubature formula involves the minimal number of points, 
which all he in the domain 12. 


Example 6.2 (Barycentric Wachspress coordinates on a pentagon). 
Here we consider the pentagon C of vertices Vi = (0,1), ^2 = (1,0), Us = 
(-l,0),n4 = (-0.5,-l),n5 = (0.5,-l). 
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To this pentagon, we associate (Wachpress) harycentric coordinates 
IHI, which are defined as follows. The weighted function associated 
to the vertex vt is defined as: 


A-i-i — Bi -\- Ai 

where is the signed area of the triangle (x, Uj) and is the 
area of (x, Vj+i, The coordinate function associated to Vi is: 


Wi(x) 


Ai(x) 


Wi(x) 


These coordinate functions satisfy: 

• Ai(x) > 0 for X E C 

• ELi a* = 1, 

• ■ Ai(x) = X 

For all polynomials p E R = M[mo, Mi, M 2 , Ms, we consider 



p o A(x)(ix 


We look for a cubature formula a E R* of the form: 


(8) I P) = 

i=i 


with Wi > 0, Ci ^ such that I[p] = {a \ p) for all polynomial p of 
degree < 2. 
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The moment matrix associated to B = {1,mq,mi,M 2 ,Ms,M 4 } 

involves moments of degree ^ 2 ; 


/ 2.5000 

0.5167 

0.5666 

0.5167 

0.4500 

0.4500 \ 

0.5167 

0.2108 

0.1165 

0.0461 

0.0440 

0.0992 

0.5666 

0.1165 

0.2427 

0.1165 

0.0454 

0.0454 

0.5167 

0.0461 

0.1165 

0.2108 

0.0992 

0.0440 

0.4500 

0.0440 

0.0454 

0.0992 

0.1701 

0.0911 

0.4500 

0.0992 

0.0454 

0.0440 

0.0911 

0.1701 ) 


Its rank is rank(if^’^) = 5. 

We compute . In this matrix there are 105 unknown param¬ 

eters. We solve the following SDP problem 


(9) 


min trace 
s.t. 0 


which yields a solution with minimal rank 5. Since the rank of the solu¬ 
tion matrix is the rank of , we do have a flat extension. Applying 
Algorithm^ we find the orthogonal polynomials pi, the matrices of the 
operators of multiplication by a variable, their common eigenvectors, 
which gives the following cubature points and weights: 


Points 

Weights 

(0.249888, -0.20028,0.249993, 0.350146, 0.350193) 
(0.376647,0.277438, -0.186609, 0.20327, 0.329016) 
(0.348358,0.379898,0.244967, -0.174627, 0.201363) 
(-0.18472,0.277593,0.376188, 0.329316, 0.201622) 
(0.242468,0.379314,0.348244,0.200593, -0.170579) 

0.485759 

0.498813 

0.509684 

0.490663 

0.51508 
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